{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>采用5折交叉验证，分别用log似然损失和正确率，对Logistic回归模型的正则超参数调优。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "#导包\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pregnants</th>\n",
       "      <th>Plasma_glucose_concentration</th>\n",
       "      <th>blood_pressure</th>\n",
       "      <th>Triceps_skin_fold_thickness</th>\n",
       "      <th>serum_insulin</th>\n",
       "      <th>BMI</th>\n",
       "      <th>Diabetes_pedigree_function</th>\n",
       "      <th>Age</th>\n",
       "      <th>Triceps_skin_fold_thickness_Missing</th>\n",
       "      <th>serum_insulin_Missing</th>\n",
       "      <th>Target</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.639947</td>\n",
       "      <td>0.866045</td>\n",
       "      <td>-0.031990</td>\n",
       "      <td>0.670643</td>\n",
       "      <td>-0.181541</td>\n",
       "      <td>0.166619</td>\n",
       "      <td>0.468492</td>\n",
       "      <td>1.425995</td>\n",
       "      <td>-0.647760</td>\n",
       "      <td>1.026390</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-0.844885</td>\n",
       "      <td>-1.205066</td>\n",
       "      <td>-0.528319</td>\n",
       "      <td>-0.012301</td>\n",
       "      <td>-0.181541</td>\n",
       "      <td>-0.852200</td>\n",
       "      <td>-0.365061</td>\n",
       "      <td>-0.190672</td>\n",
       "      <td>-0.647760</td>\n",
       "      <td>1.026390</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.233880</td>\n",
       "      <td>2.016662</td>\n",
       "      <td>-0.693761</td>\n",
       "      <td>-0.012301</td>\n",
       "      <td>-0.181541</td>\n",
       "      <td>-1.332500</td>\n",
       "      <td>0.604397</td>\n",
       "      <td>-0.105584</td>\n",
       "      <td>1.543781</td>\n",
       "      <td>1.026390</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.844885</td>\n",
       "      <td>-1.073567</td>\n",
       "      <td>-0.528319</td>\n",
       "      <td>-0.695245</td>\n",
       "      <td>-0.540642</td>\n",
       "      <td>-0.633881</td>\n",
       "      <td>-0.920763</td>\n",
       "      <td>-1.041549</td>\n",
       "      <td>-0.647760</td>\n",
       "      <td>-0.974289</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-1.141852</td>\n",
       "      <td>0.504422</td>\n",
       "      <td>-2.679076</td>\n",
       "      <td>0.670643</td>\n",
       "      <td>0.316566</td>\n",
       "      <td>1.549303</td>\n",
       "      <td>5.484909</td>\n",
       "      <td>-0.020496</td>\n",
       "      <td>-0.647760</td>\n",
       "      <td>-0.974289</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   pregnants  Plasma_glucose_concentration  blood_pressure  \\\n",
       "0   0.639947                      0.866045       -0.031990   \n",
       "1  -0.844885                     -1.205066       -0.528319   \n",
       "2   1.233880                      2.016662       -0.693761   \n",
       "3  -0.844885                     -1.073567       -0.528319   \n",
       "4  -1.141852                      0.504422       -2.679076   \n",
       "\n",
       "   Triceps_skin_fold_thickness  serum_insulin       BMI  \\\n",
       "0                     0.670643      -0.181541  0.166619   \n",
       "1                    -0.012301      -0.181541 -0.852200   \n",
       "2                    -0.012301      -0.181541 -1.332500   \n",
       "3                    -0.695245      -0.540642 -0.633881   \n",
       "4                     0.670643       0.316566  1.549303   \n",
       "\n",
       "   Diabetes_pedigree_function       Age  Triceps_skin_fold_thickness_Missing  \\\n",
       "0                    0.468492  1.425995                            -0.647760   \n",
       "1                   -0.365061 -0.190672                            -0.647760   \n",
       "2                    0.604397 -0.105584                             1.543781   \n",
       "3                   -0.920763 -1.041549                            -0.647760   \n",
       "4                    5.484909 -0.020496                            -0.647760   \n",
       "\n",
       "   serum_insulin_Missing  Target  \n",
       "0               1.026390       1  \n",
       "1               1.026390       0  \n",
       "2               1.026390       1  \n",
       "3              -0.974289       0  \n",
       "4              -0.974289       1  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train=pd.read_csv(\"FE_pima-indians-diabetes.csv\")\n",
    "train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train=train[\"Target\"]\n",
    "X_train=train.drop(['Target'],axis=1)\n",
    "feat_names=X_train.columns\n",
    "#sklearn的学习器大多之一稀疏数据输入，模型训练会快很多\n",
    "#查看一个学习器是否支持稀疏数据，可以看fit函数是否支持: X: {array-like, sparse matrix}.\n",
    "#可自行用timeit比较稠密数据和稀疏数据的训练时间\n",
    "from scipy.sparse import csr_matrix\n",
    "X_train = csr_matrix(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#GridSearchCV\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "lr=LogisticRegression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>默认参数的Logistc回归</h3>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>使用log似然损失评估</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logloss of each fold is:  [0.48753578 0.53076598 0.45433853 0.42098118 0.47917703]\n",
      "cv logloss is: 0.4745597003175165\n"
     ]
    }
   ],
   "source": [
    "# 交叉验证用于评估模型性能和进行参数调优（模型选择）\n",
    "#分类任务中交叉验证缺省是采用StratifiedKFold\n",
    "#采用5折交叉验证\n",
    "from sklearn.model_selection import cross_val_score\n",
    "loss=cross_val_score(lr,X_train,y_train,cv=5,scoring='neg_log_loss')\n",
    "print ('logloss of each fold is: ',-loss)\n",
    "print ('cv logloss is:', -loss.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>使用f1评估</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logloss of each fold is:  [0.63829787 0.58490566 0.6        0.67368421 0.60869565]\n",
      "cv logloss is: 0.6211166790836027\n"
     ]
    }
   ],
   "source": [
    "# 交叉验证用于评估模型性能和进行参数调优（模型选择）\n",
    "#分类任务中交叉验证缺省是采用StratifiedKFold\n",
    "#采用5折交叉验证\n",
    "from sklearn.model_selection import cross_val_score\n",
    "loss=cross_val_score(lr,X_train,y_train,cv=5,scoring='f1')\n",
    "print ('logloss of each fold is: ',loss)\n",
    "print ('cv logloss is:', loss.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>使用accuracy评估</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "logloss of each fold is:  [0.77922078 0.71428571 0.76623377 0.79738562 0.76470588]\n",
      "cv logloss is: 0.7643663526016468\n"
     ]
    }
   ],
   "source": [
    "# 交叉验证用于评估模型性能和进行参数调优（模型选择）\n",
    "#分类任务中交叉验证缺省是采用StratifiedKFold\n",
    "#采用5折交叉验证\n",
    "from sklearn.model_selection import cross_val_score\n",
    "loss=cross_val_score(lr,X_train,y_train,cv=5,scoring='accuracy')\n",
    "print ('logloss of each fold is: ',loss)\n",
    "print ('cv logloss is:', loss.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Logistic Regression + GridSearchCV</h3>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "logistic回归的需要调整超参数有：C（正则系数，一般在log域（取log后的值）均匀设置候选参数）和正则函数penalty（L2/L1） 目标函数为：J = C* sum(logloss(f(xi), yi)) + penalty\n",
    "\n",
    "在sklearn框架下，不同学习器的参数调整步骤相同：\n",
    "\n",
    "设置参数搜索范围\n",
    "生成学习器实例（参数设置）\n",
    "生成GridSearchCV的实例（参数设置）\n",
    "调用GridSearchCV的fit方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>使用log似然损失评估</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=5, error_score='raise',\n",
       "       estimator=LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n",
       "          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,\n",
       "          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,\n",
       "          verbose=0, warm_start=False),\n",
       "       fit_params=None, iid=True, n_jobs=4,\n",
       "       param_grid={'penalty': ['l1', 'l2'], 'C': [0.01, 0.1, 1, 10, 100, 1000]},\n",
       "       pre_dispatch='2*n_jobs', refit=True, return_train_score='warn',\n",
       "       scoring='neg_log_loss', verbose=0)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#需要调优的参数\n",
    "# 请尝试将L1正则和L2正则分开，并配合合适的优化求解算法（slover）\n",
    "penaltys=['l1','l2']\n",
    "Cs = [0.01, 0.1, 1, 10, 100, 1000]\n",
    "turns_para=dict(penalty=penaltys,C=Cs)\n",
    "lr_penalty= LogisticRegression(solver='liblinear')\n",
    "grid= GridSearchCV(lr_penalty, turns_para,cv=5, scoring='neg_log_loss',n_jobs = 4,)\n",
    "grid.fit(X_train,y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4742320712121342\n",
      "{'C': 1, 'penalty': 'l1'}\n"
     ]
    }
   ],
   "source": [
    "print(-grid.best_score_)\n",
    "print(grid.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "d:\\Anaconda3\\lib\\site-packages\\sklearn\\utils\\deprecation.py:122: FutureWarning: You are accessing a training score ('std_train_score'), which will not be available by default any more in 0.21. If you need training scores, please set return_train_score=True\n",
      "  warnings.warn(*warn_args, **warn_kwargs)\n",
      "d:\\Anaconda3\\lib\\site-packages\\sklearn\\utils\\deprecation.py:122: FutureWarning: You are accessing a training score ('mean_train_score'), which will not be available by default any more in 0.21. If you need training scores, please set return_train_score=True\n",
      "  warnings.warn(*warn_args, **warn_kwargs)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VNX5+PHPM5M9BAgQtoQQQFAIS6IRXCoWVEBUVtlkcan121+xrrWitlbRtlqX2la/XxcqWkUiLqyiuCBFqSyRfZHFsIUt7IEEss3z+2MmGELIDCSTSTLP+/WaV+4995yZ5waYh3PPPeeKqmKMMcZUxBHoAIwxxtR8liyMMcZ4ZcnCGGOMV5YsjDHGeGXJwhhjjFeWLIwxxnhlycIYY4xXliyMMcZ45ddkISL9RGSjiGwRkQlnqTNcRNaLyDoRea9UebGIrPS8ZvkzTmOMMRUTf83gFhEnsAm4DsgClgGjVHV9qTrtgWlAb1U9LCJNVTXbc+y4qtbz9fOaNGmiSUlJVXkKxhhT533//fcHVDXOW70QP8bQHdiiqpkAIpIODATWl6rzS+AVVT0MUJIozkdSUhIZGRmVCNcYY4KPiGz3pZ4/L0PFAztL7Wd5ykrrAHQQkUUislhE+pU6FiEiGZ7yQX6M0xhjjBf+7FlIOWVlr3mFAO2BnwMJwDci0llVjwCJqrpbRNoC80Vkjar+eNoHiNwF3AWQmJhY1fEbY4zx8GfPIgtoVWo/AdhdTp2ZqlqoqluBjbiTB6q62/MzE1gApJb9AFV9XVXTVDUtLs7rJTdjjDHnyZ89i2VAexFpA+wCRgK3lKkzAxgFvCUiTXBflsoUkVggT1XzPeVXAn/1Y6zGmFqssLCQrKwsTp48GehQaqyIiAgSEhIIDQ09r/Z+SxaqWiQidwPzACfwpqquE5GJQIaqzvIc6yMi64Fi4CFVPSgiVwCviYgLd+/nmdJ3URljTGlZWVnExMSQlJSESHlXwIObqnLw4EGysrJo06bNeb2HP3sWqOpcYG6ZssdLbSvwgOdVus5/gS7+jM0YU3ecPHnSEkUFRITGjRuzf//+834Pm8FtjKkTLFFUrLK/H0sWxpigNOK17xjx2neBDqPWsGQBMPkG98sYY85TvXo/LTjRr18/GjZsyI033lhu3fHjx5OSkkKnTp2IjIwkJSWFlJQUPvzww3P6zOXLl/PZZ59VKm5f+XXMwhhjgtFDDz1EXl4er732WrnHX3nlFQC2bdvGjTfeyMqVK8/rc5YvX87atWvp16+f98qVZD0LYN2eo6zbczTQYRhj6ohrrrmGmJiY82q7efNm+vbtyyWXXELPnj3ZtGkTAOnp6XTu3Jlu3brRq1cvTpw4wcSJE5kyZcp59UrOlfUsjDF1ypOz17F+d47Xeuv3uOv4Mm7RqWV9/nhTcqVj88Vdd93FpEmTaNeuHYsWLeLuu+/m888/58knn2TBggU0a9aMI0eOEBkZyeOPP87atWt56aWX/B6XJQtjjKkhjhw5wuLFixk6dOipsqKiIgCuvPJKxo0bx7BhwxgyZEi1x2bJAjjgiiFK8gMdhjGmCvjaAyjpUbz/P5f7M5xzoqo0adKk3DGMN954gyVLljBnzhy6devG6tWrqzW2oB+z2HYgl9uP/5r5hdXTxTTGmLOJjY2lRYsWTJ8+HQCXy8WqVasAyMzM5LLLLuOpp54iNjaWXbt2ERMTw7Fjx6oltqBPFkkhB/ksbAJHC0Px14OgjDHB5aqrrmLYsGF89dVXJCQkMG/ePJ/bpqen8+qrr9KtWzeSk5OZM2cOAPfffz9dunShS5cuXHvttXTu3JnevXuzatUqUlNTbYDb72JaEO84QE9dzaqso6S0ahjoiIwxtdDx48dPbX/zzTc+tUlKSmLt2rWnlbVt27bc5DJr1plPl46Li6u2h75ZsnCGcszRgGt0OX9ZtIqUkVcHOiJjTDWoSWMVtUHQX4YCOO6oT6gUE7r+Q46dLAx0OMYYU+NYsgAuiG9Gbkgsg/iamSt2BTocY4ypcSxZeETFNKKjYyfLFn8d6FCMMabGsWThIdFxFDvCuPjgJ6zJsqU/jDGmNEsWJZwhFF94IwOd/2Xa4s2BjsYY42+22vQ5sWRRSljaWBpKLrmrZ5ObXxTocIwxtUh1L1E+ffp0nnvuuUrH7Su7dRbg9k/cP13FFES3ZEDO18xe9QtGdk8MbFzGmFqpqpYoLyoqIiSk/K/pwYMHV02wPrKeRWkOJ6EX38JVzjV8/t3yQEdjjKmlKrNE+c9+9jMee+wxevbsycsvv8zMmTPp0aMHqamp9OnTh+zsbAAmTZrEfffdB8CYMWO49957ueKKK2jbtu2p5UKqkvUsypDU0Ti/eZ6Lsj9h3e5eJLdsEOiQjDHn4tMJsHeN93p7PQvx+TJu0bwLXP9M5eI6Bzk5OSxcuBCAw4cPM2DAAESEV199lRdeeIFnn332jDbZ2dksWrSINWvWMHz48CrveViyKKtRW4oSLmf4zoX8a8kOnhrcJdARGWOCzMiRI09t79ixg+HDh7N3717y8/Pp0KFDuW0GDRqEiNC1a1d27ar6+WKWLMoRkjaOpKz/x/aV88m7oSNRYfZrMqbW8LUHUNKjKBmzrEGio6NPbY8fP55HH32U/v378+WXX/LMM+WfX3h4+KltfyyKamMW5ek0kOKQaG4ons+c1XsCHY0xJogdPXqU+Ph4VJW33347YHFYsihPWDSOzoO5KWQxHy/ZFOhojDG1TGWWKC/riSeeYPDgwVx99dU0a9asCqM8N1JXnuGQlpamVbpU7/bvYHI/Hiz4Fb+85zEual6/6t7bGFOlNmzYQMeOHc+tUQ2+DOUv5f2eROR7VU3z1tZ6FmeTeBnFsW0ZEfIf0pfuDHQ0xpiqdvsnQZUoKsuSxdmI4EwdTXfHBpYuz+BEQXGgIzLGmICxZFGRbqNQcdCvaD5z19hAtzEmeFmyqEiDeGjXmxGh3zJt6dZAR2OMMQFjycILSRlNMz1A6M5FbN53LNDhGGNMQFiy8ObC/rgiGjIi5D9MtYFuY+qM2z+7nds/uz3QYdQaliy8CY3A0WUY/ZzL+GL5D5wstIFuY8yZSpYoX7lyJZdffjnJycl07dqV999//4y6VbFEOcDy5cv57LPPqiR+b2wdC1+kjiZ02RtcXbiQeet6MDAlPtARGWNqqKioKP7973/Tvn17du/ezSWXXELfvn1p2LDhqTq+LlHuzfLly1m7di39+vWrktgrYj0LX7RIQZslMzr8W95bsiPQ0RhjarAOHTrQvn17AFq2bEnTpk3Zv3+/z+03b95M3759ueSSS+jZsyebNrlXkUhPT6dz585069aNXr16ceLECSZOnMiUKVPOq1dyrqxn4QsRJGUMHec9wqFtq/lxfxfaxdXz3s4YU+2eXfosPxz6wWu9kjq+jFtc1OgiHu7+8DnHsnTpUgoKCmjXrp3Pbe666y4mTZpEu3btWLRoEXfffTeff/45Tz75JAsWLKBZs2YcOXKEyMhIHn/8cdauXctLL710zrGdK+tZ+KrrcNQRwoiQ//D+MhvoNsZUbM+ePYwdO5bJkyfjcPj2VXvkyBEWL17M0KFDSUlJYfz48ezevRuAK6+8knHjxjFp0iRcLpc/Qy+X9Sx8Fd0E6dCP4ZsXcW3GNh7s04HwEGegozLGlOFrD6CkRzG53+QqjyEnJ4cbbriBp59+mssuu8zndqpKkyZNyh3DeOONN1iyZAlz5syhW7durF69uipD9sp6FucidQz1iw/T7eRSPl+3L9DRGGNqoIKCAgYPHsy4ceMYNmzYObWNjY2lRYsWpx6L6nK5WLVqFQCZmZlcdtllPPXUU8TGxrJr1y5iYmI4dqx65n9ZsjgXF1yHRjdlbMS3TF1qA93GmDNNmzaNhQsX8tZbb526JfZc7nZKT0/n1VdfpVu3biQnJzNnzhwA7r//frp06UKXLl249tpr6dy5M71792bVqlWkpqbW7gFuEekH/B1wApNU9YxHPInIcOAJQIFVqnqLp/xW4Peeak+rauCe+lHCGYJ0G8mV373C/T9msu1AF5KaRHtvZ4yp844fPw7AmDFjGDNmjE9tkpKSWLt27Wllbdu2Lff5F7NmzTqjLC4ujip9NEMF/NazEBEn8ApwPdAJGCUincrUaQ88AlypqsnAfZ7yRsAfgR5Ad+CPIhLrr1jPSeoYnFrMkJBFpNtAtzG11uR+k/0yXlFX+fMyVHdgi6pmqmoBkA4MLFPnl8ArqnoYQFWzPeV9gS9U9ZDn2BeA/2ed+CLuQohP49bIRXyYsYOCouq/K8EYY6qbP5NFPFD6v95ZnrLSOgAdRGSRiCz2XLbytS0icpeIZIhIxrlMeqm01NEkFG6lRd5GvtxgA93GmLrPn8lCyikr+wzXEKA98HNgFDBJRBr62BZVfV1V01Q1LS4urpLhnoPOQ9GQCG6P/MYGuo0xQcGfySILaFVqPwHYXU6dmapaqKpbgY24k4cvbQMnogHS8Sb6y39Zunk3Ow/lBToiY4zxK38mi2VAexFpIyJhwEig7HD+DKAXgIg0wX1ZKhOYB/QRkVjPwHYfT1nNkTqGiOJj9HVmkL7MehfG1Dbbx45j+9hxgQ6j1vBbslDVIuBu3F/yG4BpqrpORCaKyABPtXnAQRFZD3wNPKSqB1X1EPAU7oSzDJjoKas5knpCg0TuivmOaRlZFBbbQLcxway6lyifPn06zz33XJXF741f51mo6lxgbpmyx0ttK/CA51W27ZvAm/6Mr1IcDkgZRfJ//krIyd3M/yGbvsnNAx2VMSbAqnKJ8qKiIkJCyv+aHjx4cNUHXwGbwV0ZKbcgKOOi/2sD3cYYoPJLlP/sZz/jscceo2fPnrz88svMnDmTHj16kJqaSp8+fcjOds8wmDRpEvfddx/gngh47733csUVV9C2bdtTy4VUJVtIsDJikyDpKkbt/Ya/brqRrMN5JMRGBToqY4La3j//mfwN3pcoP/mDu44v4xbhHS+i+aOPnnMs57NEObgXIly4cCEAhw8fZsCAAYgIr776Ki+88ALPPvvsGW2ys7NZtGgRa9asYfjw4VXe87CeRWWljKbhySy6yw9My8gKdDTGmBrifJYoLzFy5MhT2zt27KBPnz506dKFF198kXXr1pXbZtCgQYgIXbt2ZdeuXZWKvTzWs6isTgNg7kPcHb2Eh5alck/vCwhxWg42JlB87QGU9Chav/PvKo/hfJcoLxEd/dOac+PHj+fRRx+lf//+fPnllzzzzBlL7AEQHh5+ats9HFy17FutssKiofNgLj/5LTk5R1iwsRpnkhtjapzKLFFenqNHjxIfH4+q8vbbgVtP1ZJFVUgZQ0hxHiOjv7eBbmOCXGWXKC/riSeeYPDgwVx99dU0a9asCiM9N+KP7kogpKWlaXUt1XsGVXj5UnYWRHP1gd+xaEJvWjSIDEwsxgShDRs20LFjx3Nq48/LUDVVeb8nEfleVdO8tbWeRVUQgZRbaHVsJYnsYdoyG+g2pqZr/c6/gypRVJYli6rSbRSIg/ubZPD+sh0Uu+pGj80YY8CSRdWp3wIuuJY+RfPZezSPhZtsoNuY6lRXLqn7S2V/P5YsqlLKaCJP7KN/1A820G1MNYqIiODgwYOWMM5CVTl48CARERHn/R42z6IqXXg9RMby68jF3PRDMvtyTtKs/vn/4RhjfJOQkEBWVtY5LasRbCIiIkhISDjv9pYsqlJIOHQZTseMydRzDeeDjJ3c3bt9oKMyps4LDQ2lTZs2gQ6jTrPLUFUtdTTiKuC+ZqtIX7YTlw10G2PqAEsWVa1FN2jWhcGygKzDJ/h2y4FAR2SMMZVmycIfUsfQ8Mg6Lo3cbQPdxpg6wZKFP3QZBo5Qfts0gy/W7yP72MlAR2SMMZViycIfohvDhdeTlvMF4irkw+9tRrcxpnazZOEvqWNwnjjI/7TYzPs20G2MqeUsWfhLu2ugXnNGh3/L9oN5fJd5MNARGWPMebNk4S/OEOg2gub7FtI24jjv2UC3MaYWs2ThTyljEC3mkfjVfL5uLweP5wc6ImOMOS+WLPwprgMkdKdn3ucUFrv4aLkNdBtjaidLFv6WOprww5sY2XI/U5futIXOjDG1kiULf0seAiGR3BXzX7YeyGVx5qFAR2SMMefMkoW/RdSHTgNos/cz4iJcpC+zgW5jTO1jyaI6pIxG8nOYkLSFT9fs5XBuQaAjMsaYc2LJojokXQUNE+lX+BUFNtBtjKmFLFlUB4cDUkYTvetbrosvYOrSHTbQbYypVSxZVJduowDl3sYZ/Lg/l4zthwMdkTHG+MySRXWJbQ1tetIpezb1wx1MXWID3caY2sOSRXVKGYPjyHbuvSCbOWv2cCTPBrqNMbWDJYvq1PEmCK/PUMd/KChyMX3FrkBHZIwxPrFkUZ3CoqDzEBpu+5TLWoaSbjO6jTG1hCWL6pYyBgrzuD9+PRv3HWP5jiOBjsgYY7yyZFHdEtKgSQfSDs8lOsxpz+g2xtQKliyqmwikjMaZtYQ7Lipmzurd5JwsDHRUxhhToXNOFiLiEJH6/ggmaHQbCeJkbOS3nCx0MdMGuo0xNZxPyUJE3hOR+iISDawHNorIQ/4NrQ6LaQ4XXEvTzOl0aRHNlCU2o9sYU7P52rPopKo5wCBgLpAIjPXWSET6ichGEdkiIhPKOX6biOwXkZWe152ljhWXKp/lY5y1R+oYOLaH+9tm8cPeY6zKOhroiIwx5qx8TRahIhKKO1nMVNVCoML/CouIE3gFuB7oBIwSkU7lVH1fVVM8r0mlyk+UKh/gY5y1R4d+ENWYq3I/JzLUSboNdBtjajBfk8VrwDYgGlgoIq2BHC9tugNbVDVTVQuAdGDg+QZa54SEQZfhhG7+lJHJ0cxatZtjNtBtjKmhfEoWqvoPVY1X1f7qth3o5aVZPLCz1H6Wp6ysoSKyWkQ+FJFWpcojRCRDRBaLyCBf4qx1UkdDcQF3NMwgr6CYWat2BzoiY4wpl68D3Pd6BrhFRP4lIsuB3t6alVNW9tLVbCBJVbsCXwJvlzqWqKppwC3ASyLSrpy47vIklIz9+/f7cio1S/Mu0LwrCds+5qLmMTbnwhhTY/l6GeoOzwB3HyAOuB14xkubLKB0TyEBOO2/zqp6UFXzPbtvAJeUOrbb8zMTWACklv0AVX1dVdNUNS0uLs7HU6lhUscge1czvuMJ1u7KYY0NdBtjaiBfk0VJL6E/MFlVV1F+z6G0ZUB7EWkjImHASOC0u5pEpEWp3QHABk95rIiEe7abAFfivmW37ukyDJxh9Cn4iohQB1PtGd3GmBrI12TxvYh8jjtZzBORGMBVUQNVLQLuBubhTgLTVHWdiEwUkZK7m+4RkXUisgq4B7jNU94RyPCUfw08o6p1M1lENYIL+xO+4UMGdI5j5opd5OYXBToqY4w5jfgyGUxEHEAKkKmqR0SkMRCvqqv9HaCv0tLSNCMjI9BhnJ/NX8CUm/mx1/9xzacNeGZIF0Z2Twx0VMaYICAi33vGhyvk691QLtxjDr8XkeeBK2pSoqj12vWGmBa0zZpBh2b1mLpsp/c2xhhTjXy9G+oZ4F7c4wbrcV8++os/AwsqDid0G4ls+YI7ukWyaucR1u22gW5jTM3h65hFf+A6VX1TVd8E+gE3+C+sIJQyBtTFQPmGsBAH6Uutd2GMqTnOZdXZhqW2G1R1IEGvyQXQ6jIi16VzQ+fmzFixi7wCG+g2xtQMviaLvwArROQtEXkb+B74s//CClKpo+HAJu5se4hj+UV8snpPoCMyxhjA9wHuqcBlwMee1+Wqmu7PwIJS8mAIjaLT3lm0i4u2Gd3GmBqjwmQhIheXvIAWuGdl7wRaespMVQqPgU4DkbUfMeaSOJbvOMLGvccCHZUxxhDi5fgLFRxTvK8PZc5VymhYNZVhUSv5i7MRU5fu4IkByYGOyhgT5CpMFqrqbWVZU9VaXwkNW1NvQzr9Ok/k4+VZTLj+IiJCnYGOzBgTxHydZzGknNc1ItLU3wEGHYfD/RS9rQu5tZOQc7KIuWtsoNsYE1i+3g31C2ASMNrzegN4AFgkIl4fr2rOUbdRgHDxoU9JahxlA93GmIDzNVm4gI6qOlRVh+J+TGo+0AN42F/BBa2GraDt1ciq9xh1aQLLth1mS7YNdBtjAsfXZJGkqvtK7WcDHVT1EFDrnwV6+2e3c/tntwc6jNOljIEjOxjZdDuhTmGqzeg2xgSQr8niGxGZIyK3isituJ9LsVBEooEj/gsviHW8EcIb0OCHafTp1JyPlmdxsrA40FEZY4KUr8liPDAZ9zLlqbgffzpeVXPtjik/CY2EzkNg/SzGpMRyJK+Qeev2BjoqY0yQ8nUGtwLfAvNxPyt7ofryIAxTOaljoegEPU78h8RGNtBtjAkcX2+dHQ4sBW4GhgNLRORmfwZmgPiLIe4iHCunMLJ7KxZnHiJz//FAR2WMCUK+XoZ6DLhUVW9V1XFAd+AP/gvLACDintGdtZQRSScIcQjp9mAkY0wA+JosHKqaXWr/4Dm0NZXRdQSIk8abP+Dajs348Pss8otsoNsYU718/cL/TETmichtInIb8Akw139hmVNimkH7PrAqnVGXtuRQbgFfrN/nvZ0xxlQhXwe4HwJeB7oC3YDXVdUm41WX1NFwfB9XsYr4hpE20G2MqXY+X0pS1Y9U9QFVvV9Vp/szKFNG+74Q1QTHqimMvLQVi7YcZPvB3EBHZYwJIt6eZ3FMRHLKeR0TkZzqCjLohYS5xy42fsrw5CicNtBtjKlmFSYLVY1R1frlvGJUtX51BWlwX4pyFdJs22x6X9SUDzJ2UlDkCnRUxpggYXc01RbNkqFFCqx4l1HdW3HgeAFfbbCBbmNM9bBkUZukjoF9a7g6Zg8tGkTwng10G2OqiSWL2qTzUHCG4Vz1HiMubcW3Ww6w81BeoKMyxgQBSxa1SVQjuOhGWDON4SlNEeB9G+g2xlQDSxbAsYJj1Jp1EVNHw4nDtNy3gJ9f2JRpGTspLLaBbmOMfwV9sth6dCubDm5k7w/LWXtgbaDD8a5tL6gfDyunMKp7ItnH8pn/Q7b3dsYYUwlBnyyS6icxekkoxyKVWz65hacXP01OQQ2eQuJwQreRsOVLerUopFn9cNJtoNsY42dBnyxEhJSdITz0WQSjLhrFB5s+YMD0AXyS+UnNvTSVMhrURcjaaYxIa8WCTfvZdeREoKMyNUyPyUPpMXlooMOoVnbO/hP0yaJEZKHwSI9HeO+G92gR3YIJ30zgl1/8kq1HtwY6tDM1bgeJl8PKKQxPSwBsoNubGvmcdWNqEUsWAKV6EMmNk3m3/7s81uMx1h9Yz9BZQ3l5xcucLDoZwADLkToGDm4h4fhaeraPY9qynRTZQPdZrd+Tw/o9NfjyojE1XNAni+LjuTTblUfUscJTZU6Hk5EXjWTW4Flc1/o6Xlv9GkNmDeHbXd8GMNIyOg2C0GhY+S6juieyN+ck/9m0P9BRGWPqqKBPFq5c92NKm2SfZPfDEyg+/tNjS5tENuHZns/yRp83cIqT//fl/+PBBQ+yL7cGLLMRXg+SB8Ha6VzTLpq4mPBzWrp8xGvfMeK17/wYoDGmLgn6ZBHarBn74qM4GhvG0dmz2Tp4CCdWrjytzmUtLuOjAR8xPmU8C3YuYODMgby7/l2KXEUBitojZTQUHCN04xyGXZLA/B+y2XPUBrqNMVUv6JMFACIcbRRO63ffgeJito0ew/5XXkGLfkoGYc4wftXtV8wYOIOUpik8u+xZbvnkFlbvXx24uFtfAbFtYOUURl6aiEth2rIsn5puC3uebWHP+zlAY0xdYcmilKiLL6bNzBnUv/56DvzzZbaPu5WCrF2n1WlVvxX/d83/8fzVz3PwxEHGzB3DU989xdH8o9UfsIi7d7HtGxJlH1e1b8K0jJ0Uu2roLb/GmFrLr8lCRPqJyEYR2SIiE8o5fpuI7BeRlZ7XnaWO3Soimz2vW/0ZZ2nOmBjin3+Ols/9lfyNG9k6aBBHZ88pGzd9k/oyc9BMRncczYebP2TAjAHM/nF29c/NSBkFCKyayqjuiew6coKFm22g2xhTtfyWLETECbwCXA90AkaJSKdyqr6vqime1yRP20bAH4EeQHfgjyIS669Yy9PgpptoM3MG4e3bs/uhh9j1u99RfOzYaXXqhdXj4e4Pk35DOvH14nn020e58/M7yTyaWY2BJkC7XrByKtdeFEfj6DCmLrEZ3caYquXPnkV3YIuqZqpqAZAODPSxbV/gC1U9pKqHgS+Afn6Kk/TfJJP+m+QzysMSEmj9zr9pcvfd5Mz5hK2DBpO3fMUZ9To27sg717/DHy77AxsObWDorKH8Y/k/qm9uRspoOLqDsJ3fcnNaAl/9kE12Tg2bF2KMqdX8mSzigdLTirM8ZWUNFZHVIvKhiLQ6x7Z+JyEhxN09ntbvvgsibB8zhv3/fPm0wW9wz80YfuFwZg2axfVJ1/PGmjcYNHMQC7MW+j/Ii26EiAaw4l1GXppIsUv54HvfBrqNMcYX/kwWUk5Z2Qv6s4EkVe0KfAm8fQ5tEZG7RCRDRDL27/fvdfqoi1NpM2M6DW66kQOvvML2MWMpyDrzC7lJZBP+fNWf+VeffxHmDGP8V+N5YMED7M3d67/gQiOg882wYTZtogu5vG1jpi7dgcsGuo0xVcSfySILaFVqPwHYXbqCqh5U1XzP7hvAJb629bR/XVXTVDUtLi6uygI/G2e9erR89llaPv88+Vu2sHXgII7OmlVu3e4tuvPRTR9xT+o9LMxayMAZA/n3un/7b25G6mgoOgnrPmZUj0SyDp/g2y0H/PNZxpig489ksQxoLyJtRCQMGAmc9s0qIi1K7Q4ANni25wF9RCTWM7Ddx1NWIzS48QbazJhB+EUXsft3D7Prtw+dMfgNEOoM5Zddf8n0gdO5uNnFPJfxHCPnjGTV/lVVH1TLi6FpJ1gxhb7JzYiNCiV9mQ10G2Oqht+ShaoWAXfj/pLfAExT1XUiMlFEBniq3SMi60RkFXAPcJun7SHgKdwJZxkw0VNWY4QlxNP67bdocs9vyPn0U7YOHEQXGf+HAAAVVklEQVTe99+XW7dVTCv+95r/5cWfv8jh/MOMnTuWJ797smrnZpTMudiVQfjhLdx8SQKfr9vH/mP53tsaY4wXfp1noapzVbWDqrZT1T95yh5X1Vme7UdUNVlVu6lqL1X9oVTbN1X1As9rsj/jnNxvMpP7nftHSEgIcb/+NUlT3gWnk+1jx7H/H/84Y/Ab3HMzrmt9HbMGzWJsp7FM3zydATMGMOvHWVU3N6PrCHCEwIp3GXFpIkUu5UMb6DbGVAGbwV0FIlNSaDP9YxoMGMCB//0/to8eQ8GO8i8BRYdG89ClD/H+je+TEJPAY98+xh3z7iDzSBXMzagXB+37wqp0LmgcTvc2jUhfZgPdxpjKs2RRRZz16tHymb8Q/+IL5GdmsnXQYI7MmHHWXsOFjS7knevf4Y+X/5FNhzcxdPZQ/r7875woquRCgKmjITcbtnzJLd0T2X4wj8WZByv3nsaYoGfJoorV79+ftjNnEN6pI3smPMLuBx+kOKf8h+44xMHNHW5m1qBZ9G/Tn0lrJjF45uDKzc1o3wei42DFu/Tr3JwGkaG8Z8/oNsZUkiULPwht2ZLWb79N3H33kTPvczIHDSJv2bKz1m8c2Zg//exPvNn3TcKd4Yz/ajz3fX3f+c3NcIa6xy42fUZEwWGGXBzPvHV7OXjcBrqNMefPkoWfiNNJk1/9D0lT30NCQtl+621kv/QSWlh41jaXNr+UD2/6kHsvvpdFuxYxYMYA3l73NoWus7cpV+oYcBXB6mmM6p5IYbHy0XIb6DbGnD9LFn4W2bUrbT7+mAaDBnHw1dfYNnoMBdu3n7V+qDOUO7vcyYxBM+jevDvPZzzPiDkjWJm98qxtztC0o3vexYp36dC0HmmtY0lfurP6V8Q1xtQZliyqgbNeNC3//CfiX/obBdu2sXXwEI58PL3CL+/4evH8s/c/eanXS+Tk5zD207E88d8nOHLyiG8fmjoastfBnpWM6p5I5oFclmytUVNVjDG1iCWLalS/Xz/azpxBRHIyex59lF33P0Dx0bNPzBMRrkm8hlmDZnFb8m3M2DKDATMGMGPL2e+yOqXzUHCGw4op9O/SgpiIkHN6RrcxxpRmyaKahbZoQeJbk4l74AGOffklmQMHkbt0aYVtokKjeDDtQd6/8X1a12/NHxb9gds+u40th7ecvVFkLHS8CdZ8QKQUMiQ1nk/X7uVwbkEVn5ExJhhYsggAcTppctcvSZr6Ho7wcHbcehvZL/6twsFvcM/NePv6t3nyiif58eiPDJs9jL99/zfyCvPKb5A6Gk4egY1zGdUjkYIiFx+v2FV+XWOMqYAliwCK7NKFNh9/RIOhQzj4+utsG3ULBdu2VdjGIQ6GtB/C7EGzubHdjby59k0GzxzMgp0Lzqzc5mqonwArp3BR8/qkJjZk6tIdNtBtjDlnliwCzBEdTcunnyb+73+nYOdOMocM5chHH3n9Qo+NiOWpK5/irX5vERUaxW/m/4Z75t/DnuN7Sr250/2M7h/nw9FdjLo0kS3Zx8nYftjPZ2WMqWssWdQQ9fv2oe3MGUR26cKex37Prnvvo/iI9zufLml2CdNumsb9l9zP4j2LGThzIJPXTv5pbkbKLaAuWDWVG7u1oF64DXQbY86dJYsaJLR5cxLf/BdNf/sgx77+2j34vXiJ93aOUO7ofAczBs6gR4sevPj9iwyfPZwV2SugUVtofSWsnEJUqJNBqS35ZPUeXMVh1XBGxpi6wpJFDSNOJ43vvJOkqVNxREay4/bbyX7hBbTA+11MLeu15J+9/8k/ev2D3MJcxn06jscXPc7hLkPgUCbsWMzISxPJL3JxMqddNZyNMaausGRRQ0V2TqbNxx/RcNgwDr4xiW2jbiE/c6tPbXsl9mLGwBnc3vl2Zv84mwGb32J6w0a4VrxD5/gGdE1owIkjF2Lj3MYYX4UEOgBzdo6oKFpMfJLoq37G3t//ga1Dh9LskQk0HDYMEamwbVRoFA9c8gA3tb2Jpxc/zeMFR5mePZ/fZ69iVPdEHvn4KAcyR9D7hQVEh4UQHe4kOiyEqPAQ6oU7iQoLITrMSXS4uyw6zF1WLzyEKE/dn9o4CQ9xVtNvxRgTCJYsaoH6111HZNdu7J7wMHsf/yO533xD84kTCYmN9dq2fWx7JvebzMwlL/Li+jcZ8ek4RnUcS1SjPFzFMXRscQG5+UXk5RezN+ckeQXFHM8vIi+/iNyCYp9jDHVKuQkmuiTRnLbvTjD1wkNOaxNdkqQ8dUOc1vE1pqaQunLPfVpammZkZAQ6DL9Sl4tDb71N9t/+RkhsLC2ffYboyy/3sbFy5OWL+VtMOB9LLiEKoQodm1+MQxyICA7cPwXBIQ5PM8Gl4HKBSwWXC4o9+0UuxeUSij3bxcWc2i4q2S6GwmI9VVZUrID89NKSHpIDPbXtPuYUB6FOB2EhTsIcTkJDnISV7DvdvZmftkNO7UeEuI+Fh4QQEerefjnjbcRRzMOX/4piLcalisulFLtc7m11Uawl2yXlPx07vfz0Mve2y/170rLtQEuVKUqxKqrqKVdcuPdd6jr1s6Suej5D8RyjVFs4te2pfVq77NwDADSJalT+XwnO9m9fKX1Iyx4rt7z8ktOOlHP47DH81EDPLDrrZ+XkHweBmLDoCt63bskpyEWKo1n36w/Pq72IfK+qaV7rWbKofU6uX8+u3z5EwdatNLrjdpreey8S5sPdTQufh/lPsWLMe/x64QRcAsktL/3pi6rky8nzxXTaNqfXKfdYmXoudYGCC/eXJXDqS/mML8eS98WFKni+FqnoC6imcSc7T8IrvX0qMZbaBlCH5+x8aec4VSal6pdsi7gTr/vypPtVVFyEAuFOL383Sl3SPLWlpXfOaFBBacWXR8v5JC9NKnpXOWMv17OaQXRoVIVx1AwV/658lVuYB8XRrPv1R+cXhY/Jwi5D1UIRnTrR5qMP2ffssxz615vkfvcd8c8/T3jbthU37DYKvv4TqTtW0KrIPcbwZt83qyHiyjlbQipJQPlFReQVFHE8v5DcgkLy8os4XlBIXkExeQUF5OYX8c9lU1GXk7FdBuJwOHCK++VwCCEiOJyeMocDp7h7Vk6H4BSnu8whhDgcOEQIESdOpwOHQIjDSYjD3SNzOgSnw70ApNOz7xDB4cD9np59p6fMUbqO4N52eI6Xbudpe656TB4KwJLbz+9LpDYK5nP2N0sWtZQjMpIWTzxBvauuYs9jv2frkKE0mzCBhiOGn33wu0E8tOsNK6ci0SGol0HymqLspbGyokIhNrLi93hz858AmNDHx8t2xpjT2AhiLRdzzTW0mTWTqIsvZu8TT5B1928oOlzBch4poyEni875J6svSGNMrWfJog4IbdqUVpPeoOmEh8lduJCtAwZyfNGi8itf2B8iGvKLH45y2zRbrtwY4xtLFnWEOBw0vu02kqa9j6N+fXb+4k72PfMsrrIzv0MjoMswmjU5QbPGJ2DvGjiZE5igjTG1ht0NVQe5Tpwg+7nnOPzeVMI7diT++ecIb1dqeY996yn65xWEhJT6s4+MhYatIbZ1qZ9J7p8NWrmTTC0299pOAPT/cn2AI6k+ds7BobLnbHdDBTFHZCTNH3+c6KuuYs+jj3kGvx+m4ciR7sHvZp2Yv6w50ZFFXPmn5+Dwdjiy3f1z3zrY+CkUl+mRxLQoJ5l4ftaPdy+HboypsyxZ1GExvXoROWsmux95lL1PTuT4wm9o8aenCWnUiKJiB0ePh0Hy4DMbulxwfO/pSaTk5/b/wpoP3Muel3CEQIOEMkkk6af96LjT7uU3xtQ+lizquJC4OFq9/hqH332X7OeeJ3PgQFr+5ZmKGzkcUL+l+9W6nFtNiwvh6M7yk8nGTyF3/+n1Q6OgYeLZeyYRDaruhI0xfmHJIgiIw0GjceOI6tGD3b/9LTvvvJNG9eBkOOT+9784oqKQqCgcUdE4oqNwREcjoaFnn6/hDHU/J6PRWSYBFuTCkR3lJ5Md30F+mQH1iIblJJEk98+GibV+vMSYusCSRRCJuPBCkj74gOznnocpU2hwHHbc8YvyK4eE4IiKOvsrOupUknFGR3uSTckr2lOnDY6WyTgu8NSNiHAvcHDisDuZlE0k+3+AzZ9DUZk5IPWa/5REGiaWGS9JAKf9NTbG3+xfWZBxRETQ/A+/Z9VXU3C64IoX38GVl+d+5eb9tH3qles55v5ZmL0PLVPP5wdjOBzlJh6JjsIR1QpH1IU4IqNwRCoO8nFIHo7iYzhOHMaRcxDHpu9w5M/E4SzCEaI4QhQJdSAN48+8g6tkv14zGy8xpgpYsghSxSFCMRCV5vWOuQqpKnry5OkJJje3TOLJPbWt5SSm4gMHKczbeXoCKj7b8uiNzyiRMBeO0K04nFtwOIpwhLhwhLqTiSPUgSO6Hh1D8yhWOHD/wNItPT+k4rKzbZ/lmCCl1ojz8TNK2pyR2CqKR0r98Hx2qWNJDncP7dCTdxEsTp3zxP8JcCTVJ8lxksJq+A+RJYsglS9eFlPykYggkZE4IiOh8Zlf5OdDVdGCglJJJfdUkinO/Wm73N7Q8WO4cg5RfDyHwtzjuA7l4dpTgOtkBKiwf9umKomxNhDCAdi345sAR1J9Tp3zewsDHEn1EcKJbljo98+xZGFqHBFBwsNxhIeDDw948sXcazuBKtd/tvr0A+VcQiv3olp5l9oqW1ZSriXPbShVVnq75DZlLbVku6vUrcul65VqP3/0taDQe8rn5Z1RnTR/dB8Aek+ZF+BIqs/80X1BQrjAz59jycIEDxEkNNR7tWoIpToUqns1H2ezxABHUn0KKTnn1gGOpPoU4qiWx77Y2lDGGGO8smRhjDHGK7sMFaSeuSUJgHIW+zDGmDP4tWchIv1EZKOIbBGRCRXUu1lEVETSPPtJInJCRFZ6Xq/6M05jjDEV81vPQkScwCvAdUAWsExEZqnq+jL1YoB7gCVl3uJHVU3xV3zGGGN858/LUN2BLaqaCSAi6cBAoOyi608BfwV+68dYTJB7cnRHAPoHOI7qZOccHKrrnP15GSoe2FlqP8tTdoqIpAKtVHVOOe3biMgKEfmPiFzlxziNMcZ44c+eRXm3q5+6G1hEHMDfgNvKqbcHSFTVgyJyCTBDRJJV9bTlSkXkLuAugMTE4LmX3Bhjqps/exZZQKtS+wnA7lL7MUBnYIGIbAMuA2aJSJqq5qvqQQBV/R74EehQ9gNU9XVVTVPVtLi4OD+dhjHGGH8mi2VAexFpIyJhwEhgVslBVT2qqk1UNUlVk4DFwABVzRCROM8AOSLSFmgPZPoxVmOMMRXw22UoVS0SkbuBeYATeFNV14nIRCBDVWdV0LwnMFFEioBi4FeqeshfsRpjjKmYXyflqepcYG6ZssfPUvfnpbY/Aj7yZ2zBrlOL+oEOwRhTi9hyH8YYY7yy5T6C1OR+kwMdgjGmFrGehTHGGK8sWRhjjPHKkoUxxhivLFkYY4zxyga4TVCwW4WNqRzR8h4mXwulpaVpRkZGoMMwxphaRUS+V9U0b/XsMpQxxhivLFkYY4zxypKFMcYYryxZGGOM8cqShTHGGK8sWRhjjPHKkoUxxhivLFkYY4zxypKFMcYYr+rMDG4R2Q9sr8RbNAEOVFE4tUWwnXOwnS/YOQeLypxza1WN81apziSLyhKRDF+mvNclwXbOwXa+YOccLKrjnO0ylDHGGK8sWRhjjPHKksVPXg90AAEQbOccbOcLds7Bwu/nbGMWxhhjvLKehTHGGK8sWXiIyHMi8oOIrBaR6SLSMNAx+ZuIDBORdSLiEpE6ffeIiPQTkY0iskVEJgQ6Hn8TkTdFJFtE1gY6luoiIq1E5GsR2eD5e31voGPyNxGJEJGlIrLKc85P+uuzLFn85Augs6p2BTYBjwQ4nuqwFhgCLAx0IP4kIk7gFeB6oBMwSkQ6BTYqv3sL6BfoIKpZEfCgqnYELgPGB8Gfcz7QW1W7ASlAPxG5zB8fZMnCQ1U/V9Uiz+5iICGQ8VQHVd2gqhsDHUc16A5sUdVMVS0A0oGBAY7Jr1R1IXAo0HFUJ1Xdo6rLPdvHgA1AfGCj8i91O+7ZDfW8/DIQbcmifHcAnwY6CFNl4oGdpfazqONfIsFORJKAVGBJYCPxPxFxishKIBv4QlX9cs4h/njTmkpEvgSal3PoMVWd6anzGO7u7JTqjM1ffDnnICDllNltgHWUiNQDPgLuU9WcQMfjb6paDKR4xlmni0hnVa3ysaqgShaqem1Fx0XkVuBG4BqtI/cUezvnIJEFtCq1nwDsDlAsxo9EJBR3opiiqh8HOp7qpKpHRGQB7rGqKk8WdhnKQ0T6AQ8DA1Q1L9DxmCq1DGgvIm1EJAwYCcwKcEymiomIAP8CNqjqi4GOpzqISFzJnZsiEglcC/zgj8+yZPGTl4EY4AsRWSkirwY6IH8TkcEikgVcDnwiIvMCHZM/eG5cuBuYh3vQc5qqrgtsVP4lIlOB74ALRSRLRH4R6JiqwZXAWKC359/wShHpH+ig/KwF8LWIrMb9n6IvVHWOPz7IZnAbY4zxynoWxhhjvLJkYYwxxitLFsYYY7yyZGGMMcYrSxbGGGO8smRhzDkQkePea1XY/kMRaevZricir4nIj54VQxeKSA8RCfNsB9WkWVOzWbIwppqISDLgVNVMT9Ek3Iv9tVfVZOA2oIlnscOvgBEBCdSYcliyMOY8iNtzIrJWRNaIyAhPuUNE/tfTU5gjInNF5GZPs9FAyRpk7YAewO9V1QXgWRX3E0/dGZ76xtQI1s015vwMwf38gG5AE2CZiCzEPYs4CegCNMU9Y/xNT5srgame7WRgpWcRuPKsBS71S+TGnAfrWRhzfn4GTFXVYlXdB/wH95f7z4APVNWlqnuBr0u1aQHs9+XNPUmkQERiqjhuY86LJQtjzk95y55XVA5wAojwbK8DuolIRf8Gw4GT5xGbMVXOkoUx52chMMLz4Jk4oCewFPgWGOoZu2gG/LxUmw3ABQCq+iOQATzpWS0VEWkvIgM9242B/apaWF0nZExFLFkYc36mA6uBVcB84Heey04f4X5+xlrgNdxPajvqafMJpyePO3E/mGqLiKwB3uCn52z0Aub69xSM8Z2tOmtMFROReqp63NM7WApcqap7Pc8b+Nqzf7aB7ZL3+Bh4JEiekW5qAbsbypiqN8fzQJow4ClPjwNVPSEif8T9/O8dZ2vseUDTDEsUpiaxnoUxxhivbMzCGGOMV5YsjDHGeGXJwhhjjFeWLIwxxnhlycIYY4xXliyMMcZ49f8Boeyz/GPWomIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "test_means=grid.cv_results_['mean_test_score']\n",
    "test_stds=grid.cv_results_['std_test_score']\n",
    "train_stds=grid.cv_results_['std_train_score']\n",
    "train_means=grid.cv_results_['mean_train_score']\n",
    "# plot results\n",
    "n_Cs = len(Cs)\n",
    "number_penaltys = len(penaltys)\n",
    "test_scores = np.array(test_means).reshape(n_Cs,number_penaltys)\n",
    "train_scores = np.array(train_means).reshape(n_Cs,number_penaltys)\n",
    "test_stds = np.array(test_stds).reshape(n_Cs,number_penaltys)\n",
    "train_stds = np.array(train_stds).reshape(n_Cs,number_penaltys)\n",
    "\n",
    "x_axis = np.log10(Cs)\n",
    "for i, value in enumerate(penaltys):\n",
    "    #pyplot.plot(log(Cs), test_scores[i], label= 'penalty:'   + str(value))\n",
    "    plt.errorbar(x_axis, -test_scores[:,i], yerr=test_stds[:,i] ,label = penaltys[i] +' Test')\n",
    "    plt.errorbar(x_axis, -train_scores[:,i], yerr=train_stds[:,i] ,label = penaltys[i] +' Train')\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel( 'log(C)' )                                                                                                      \n",
    "plt.ylabel( 'logloss' )\n",
    "plt.savefig('LogisticGridSearchCV_C.png' )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>使用f1评估</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=5, error_score='raise',\n",
       "       estimator=LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n",
       "          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,\n",
       "          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,\n",
       "          verbose=0, warm_start=False),\n",
       "       fit_params=None, iid=True, n_jobs=4,\n",
       "       param_grid={'penalty': ['l1', 'l2'], 'C': [0.01, 0.1, 1, 10, 100, 1000]},\n",
       "       pre_dispatch='2*n_jobs', refit=True, return_train_score='warn',\n",
       "       scoring='f1', verbose=0)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#需要调优的参数\n",
    "# 请尝试将L1正则和L2正则分开，并配合合适的优化求解算法（slover）\n",
    "# solver：逻辑回归损失函数的优化方法，有四种算法供选择\n",
    "# ‘newton-cg’：坐标轴下降法来迭代优化损失函数 \n",
    "# ‘lbfgs’：, ‘liblinear’：牛顿法变种 \n",
    "# ‘sag’：随机梯度下降\n",
    "\n",
    "penaltys=['l1','l2']\n",
    "Cs = [0.01, 0.1, 1, 10, 100, 1000]\n",
    "turns_para=dict(penalty=penaltys,C=Cs)\n",
    "lr_penalty= LogisticRegression(solver='liblinear')\n",
    "grid= GridSearchCV(lr_penalty, turns_para,cv=5,scoring='f1',n_jobs = 4,)\n",
    "grid.fit(X_train,y_train)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6493683608115893\n",
      "{'C': 0.01, 'penalty': 'l2'}\n"
     ]
    }
   ],
   "source": [
    "print(grid.best_score_)\n",
    "print(grid.best_params_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>使用accuracy评估</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=5, error_score='raise',\n",
       "       estimator=LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n",
       "          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,\n",
       "          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,\n",
       "          verbose=0, warm_start=False),\n",
       "       fit_params=None, iid=True, n_jobs=4,\n",
       "       param_grid={'penalty': ['l1', 'l2'], 'C': [0.01, 0.1, 1, 10, 100, 1000]},\n",
       "       pre_dispatch='2*n_jobs', refit=True, return_train_score='warn',\n",
       "       scoring='accuracy', verbose=0)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#需要调优的参数\n",
    "# 请尝试将L1正则和L2正则分开，并配合合适的优化求解算法（slover）\n",
    "# solver：逻辑回归损失函数的优化方法，有四种算法供选择\n",
    "# ‘newton-cg’：坐标轴下降法来迭代优化损失函数 \n",
    "# ‘lbfgs’：, ‘liblinear’：牛顿法变种 \n",
    "# ‘sag’：随机梯度下降\n",
    "\n",
    "penaltys=['l1','l2']\n",
    "Cs = [0.01, 0.1, 1, 10, 100, 1000]\n",
    "turns_para=dict(penalty=penaltys,C=Cs)\n",
    "lr_penalty= LogisticRegression(solver='liblinear')\n",
    "grid= GridSearchCV(lr_penalty, turns_para,cv=5,scoring='accuracy',n_jobs = 4,)\n",
    "grid.fit(X_train,y_train)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.76953125\n",
      "{'C': 0.01, 'penalty': 'l2'}\n"
     ]
    }
   ],
   "source": [
    "print(grid.best_score_)\n",
    "print(grid.best_params_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以对比发现，使用Logistic回归和GridSearchCV实现的效果要优于使用默认参数。\n",
    "使用loss似然损失得到的最佳参数是  l1正则项回归系数C=1\n",
    "使用f1和accuracy得到的最佳参数是  L2正则项回归系数C=0.01"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
